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ABSTRACT 

Using analytic calculations and N-body simulations we show that in constant density 
(harmonic) cores, sinking satellites undergo an initial phase of very rapid (super- 
Chandrasekhar) dynamical friction, after which they experience no dynamical friction 
at all. For density profiles with a central power law profile, p oc r~", the infalling 
satellite heats the background and causes a to decrease. For a < 0.5 initially, the 
satellite generates a small central constant density core and stalls as in the a — case. 

We discuss some astrophysical applications of our results to decaying satellite 
orbits, galactic bars and mergers of supermassive black hole binaries. In a companion 
paper we show that a central constant density core can provide a natural solution to 
the timing problem for Fornax's globular clusters. 
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1 INTRODUCTION 

In a seminal paper, IChandrasekhaJ lll943h showed that a 
, massive particle moving through an infinite, homogeneous 
' and isotropic background of lighter particles experiences a 

force of dynamical friction given by; 

M,^ = ~A7TGHd^j^ln(^) r\l{v')dv' (1) 

at Ivl"^ VOmin/ Jg 

where Mc and v are the mass and velocity of the in-falling 
particle, M(v')dv' is the mass density of background objects 
with speeds u' ^ u' + dv' , and 6max and 6min are the maxi- 
mum and minimum impact parameters for the encounters^ . 
From here on, we refer to the massive in-falling object as a 
'Globular Cluster' (GC) and the background of lighter par- 
ticles as simply 'particles'. However, we could equally refer 
to, for example, a bar moving in a background of stars and 
dark matter. 

While equation is only strictly valid for an infinite, 
homogeneous and isotropic background, it has been shown 

* Email: justinCSphysik.unizh.ch 

^ Note that bmin ~> c an be achieved (see e.g. IWhitelll976t 
iBinnev fc TremainjligSTl p. 423), while ^max — ^ cannot. This 
is because the derivation of equation Q assumes an infinite back- 
ground; ftmax defines a scale on which the infinite background 
should be truncated. This is often, reasonably, take to be the 
radius at which the mean density falls by a factor of two or so. 



to work remarkably well for satellites orbiting in spherical 
gala xies with more general background distributions'^ (see 
e.g. IWhite"l983', 'Zaritskv fc Whitd llQSd. ICora et aLllToOTl 
and JBontckoe & van Albada ■1987^ . Such successes make 
equation of great practical value. But they beg the ques- 
tion: why has it been so successful, even when it is used so 
far beyond its expected regime of validity? Are we missing 
important physical insight into the dynamical friction pro- 
cess in spherical systems? Does Chandrasekhar fail to work 
well in some situations? 

In order to ad dress some of these issues, 
Tremaine fc Weinberd ((1981) (hereafter TW84) and 
Weinber^ ' i^S^ (hereafter W86) formulated a perturbative 
theory of dynamical friction which could be applied to 
spherical systems. Notice from equation that most 
of the dynamical friction originates from particles with 
large impact parameters: it is the accumulation of many 
long range small interactions which leads to most of the 
dynamical friction; not the large angle scattering of close 
encounters. This is why perturbative methods can be used. 
TW84 and W86 consider a general, small, perturbation to a 
single background particle; and then sum over all particles 



^ When corrected for velocity anisotropies, it has also been shown 
to wo rk well in aspherical system s (see e.g. .Binncy 1977.. StatleJ 
I199J and iPefiarrubia et alJl200d) . 
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in the system to obtain to total torque induced on the 
perturber. 

In section |5] we will briefly summarise the essence of 
this perturbation method. For now, it is important to note 
the key assumptions in the method; and the key results. 
The two main assumptions are: (i) that the perturbation is 
small and (ii) that the frequency of the perturber changes 
with time, Q.s = ^sit), faster than the perturbation can grow 
non-Unear. For most potentials of interest this second as- 
sumption is satistifed. The perturber (the GC) will lose an- 
gular momentum to the background particles as a result of 
the dynamical friction, and Slj, will then increase as the GC 
falls inwards. 

Under the above assumptions, the perturbation method 
gives us new physical insight into the dynamical friction 
problem. To the order of the perturbation approximation, 
all of the torque comes from background particles which are 
close to resonance with the perturber. Non-resonant parti- 
cles do not contribute to the friction at all. This is a key dif- 
ference between the perturbation solution and that of equa- 
tion It suggests that if equation is ever going to fail, it 
would do so for background particle distributions which are 
especially resonant. 

In this paper, we describe such a super-resonant po- 
tential: that of the constant density (harmonic) core. For 
this special potential, all particles and the perturber always 
move with constant angular frequency, fl. In this case, per- 
turbation methods can no longer be used. This is because 
Qs = const; assumption (ii), above, is violated; and the per- 
turbations, which are always driven at the same resonant 
frequency, can grow indefinitely^. 

To cope with this special case, we develop a non- 
perturbative analytic model using a 3D driven harmonic os- 
cill ator. T his essentially generalises an earlier result derived 
bv iKalna is (1972). Using our analytic model and N-body 
simulations, we show that in constant density cores, equa- 
tion fails. Sinking satellites undergo an initial phase of 
very rapid (super-Chandrasekhar) dynamical friction, after 
which they experience little or no dynamical friction at all. 

Weinberg & Katz (2005) and Weinberg & Katz (2006), 
find similar stalling results for galactic bars (which may be 
thought of as two diametrically opposed satellites) inside 
constant density cores. 

Constant density cores have recently become interest- 
ing in astrophysics. Observations of galaxies on all scales 
from dwarf spheroidals in the Local Group, up to gi- 
ant spirals suggest that their central dark matter den- 
sity has such a constant density core on the scale of 
~ Ikpc (see e.g. 'Klcvna ct al." 2003, "do Blok et alj 1200 ll 
iBorriello fc Salucci 2001 and B in.nev & Evans 2001); but see 
also lHavashi et all EoO^ and lRhee e t al. (2004), for a dis- 
cussion of the potential systematic errors in such observa- 
tions. If cores are present at the centre of galaxies, their 
resonant properties can significantly affect the dynamics. 
Bars can be much longer lived'', while in- falling satellites 
and GCs will stall at the core radius. In a companion paper. 



^ This is true for any pert ur bative scheme fejgJCol£^^aL^^99|K 
^ iDebattista fc Sellwoodl <1998D and lDebattista~ SellwoodI 
hood) show that low central dark matter densities lead to bars 
which remain fast. Here we discuss the extreme case of constant 
density cores, in which we show that bars would not slow down 



iGoerdt et al.l (l2 006|l . we investiga te t his last idea further (see 
also [ Hernandez fc GilmorelllQQsl and lSanchez-Salcedo et all 
I2OO6I) . The Fornax dwarf spheroidal galaxy in the Local 
Group has 5 GCs at a range of projected radii. Application 
of Chandrasekhar dynamical friction suggests the clusters 
should rapidly fall to the centre of Fornax from their cur- 
rent positions; fine tuning is required to have t hem arrive at 
their present positions at the current epoch. In lGoerdt et alJ 
(2006) we show that a small core of radius greater than 
0.24 kpc can solve this problem by causing some, or all, of 
Fornax's GC to stall. 

This paper is organised as follows: in section|5|we briefly 
review the perturbative method for calculating dynamical 
friction and demonstrate that it fails for the special case of 
a constant density core. We show that insight can be gained 
from a non-perturbative approach by modelling the system 
as a driven harmonic oscillator. In sectionj^we describe our 
semi-analytic and full N-body simulations. In section 2] we 
test our analytic model against these high resolution (~ 10^ 
particles) simulations of satellites sinking in harmonic cores. 
We demonstrate that such high resolution is required in or- 
der to reduce numerical precession of the GC orbit plane, 
but that near-converged results for the GC orbit can be ob- 
tained at lower resolution with 0(1Q® particles). We discuss 
the importance of the initial GC orbit, mass, the underlying 
gravitational potential and the particle-particle interactions. 
Finally, in section |K| we briefly discuss the implications of 
these results and present our conclusions. 



2 ANALYTIC RESULTS 

2.1 A brief review of the perturbation method 

The essence of the TW84 perturbative approach to dy- 
namical friction can be understood in the following way: 
consider a spherical potential, ^{r), to which a small non- 
axisymmetric perturbation, <I>s, is applied. The perturbation 
rotates with angular frequency Q3. In this case, the equa- 
tions of motion of a test particle moving in a frame station- 
ary with respect to the perturbation are given by: 



r + V[$ + ^s] + 2n^ X r + n^ X (£2^ x r) 







(2) 



where the third and fourth terms are the familiar Coriolis 
and Centrifugal inertial forces respectively. 

The problem is symmetric about the plane containing 
the perturbation, so it makes sense to work in cylindrical 
coordinates: r = r{R, 4>, z). Equation |5| then reduces to: 



uR 
1 (9$ 

R(l> + 2R(l>= --—^ -2Rns 
R dtp 

and equation 01 can be rearranged to give: 



(3) 
(4) 

(5) 



at all. This agrees well with earlier f indinRS bv lWeinberg fc Katd 
hood) and [Weinberg fc Katd <2006l) . 
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where Jz is the 2-component of the specific angular momen- 
tum of the test particle and we have introduced the notation 
fis ~ IHsl; and similarly for other vectors. 

In order to solve equation |S1 we must now specify the 
perturbation, $s(7?, and the angular motion of the test 
particle, (pit). While it is not necessary in general, it also 
greatly simplifies the analysis to assume that _R = 0, which 
we do from here on. With this assumption, we can still il- 
lustrate usefully the key points of the perturbation method. 

We consider the perturbation: "l>s — yle"""^. This is in- 
structive since it is then one component of a more general 
Fourier series sum. We can find (j){t) if we assume that the 
perturbation is small. The usual trick is to suppose that over 
short times the particle trajectory is the same as in the un- 
perturbed case. For the unperturbed case, "I>s,r2s ^0 and 
equation |5] gives: (jtin ~ ^*t + const; where the subscript 
in reminds us that this is now with respect to an inertial 
frame. Transforming (j>in to the non-inertial frame rotating 
with Qs, gives: (p = (Q, — Q,s)t + const. 

Equation |S] may now be integrated to give: 

It is clear from equation |H] that Jz just oscillates with 
no time averaged change^ (i.e. no dynamical friction) unless 
SI, — ^la. At this resonant frequency the test particle ap- 
pears to have a pathological specific angular momentum. In 
practise this just means that the approximation that the 
perturbation is small fails. 

TW84 show that if Qs = f2s(t), then this problem 
can be solved. Provided fis changes faster than the time 
taken for the perturbation to grow into the non-linear 
regime, then we can sum over all of the resonant interac- 
tions from the background particles and calculate the re- 
sulting torque on the perturber®. There are two regimes of 
interest: fast and slow passages through resonance. The fast 
pass age through resonance recovers the LBK torque formula 
ijLvndcn-BcU & Kalnajs 197^). This is the perturbation the- 
ory equivalent of equation^ it describes the dynamical fric- 
tion. For slow passages through resonance, TW84 find quite 
different behaviour. The torque is stronger than in the LBK 
case, reversible, and can lead to the capture (gravitational 
binding) of background particles by the perturber. These 
differences led TW84 to refer to this as dynamical feedback, 
rather than friction. We return to this effect in section 

In this paper we discuss a special potential of interest 
generated by a constant density core. For this special case, 
the divergence in equation |H| persists because fls stays fixed. 
The potential for a constant density core is the harmonic 
potential given by: 

$ = —r^ + const (7) 

where fl is the angular frequency of test particles (including 
the GC; ila = f2) in the harmonic core. 

The equation of motion for the GC perturber moving 
in the harmonic potential is given by: 

^ Recall that we have assumed that _R = const. 

^ Note that there is now an extra term which should also be 

included in equation |^ S7 X r; we assume that this is small. 



f c + V$ = = ii, + f^V^ (8) 

Which may be trivially solved to give the general solution: 

= [Xsin(f7t + <;/i^),ysin(f7f + <^y)] (9) 

From equation |5] we can see that orbits in harmonic 
potentials are of fixed relative phase angle, closed and of 
constant angular frequency, SI. This means that provided 
the potential remains harmonic, any perturbation to the GC 
orbit - including dynamical friction and loss of angular mo- 
mentum - will not change S7 or Sis = S7. In other words, 
S7s 7^ f2s (t) and we can no longer apply perturbation theory 
methods. 

M. Weinberg (private communication) has made the 
valid point that the perturber itself, and the non-spherically 
symmetric background distribution it induces, cause devia- 
tions from true harmony. It may be possible to use a per- 
turbative approach in this, more realistic, case. 

2.2 A non-perturbative approach 

Perturbation methods fail for the harmonic core. However, 
all is not lost analytically. We can still gain much insight 
by writing down the equations of motion for the GC and 
a tracer background population, and searching for stable 
solutions. As we shall show next, for the special case of a 
harmonic potential plus point mass perturber (the GC), so- 
lutions exist where the background particles rotate about 
the GC on stable epicycles. Stable orbits mean no time av- 
eraged angular momentum transfer and, therefore, no dy- 
namical friction. 

Such a model allows us to make firm qualitative (if not 
quantitative) statements about what will happen when a 
GC is introduced to an isotropic constant density core. Ini- 
tially, particles will be in equilibrium in the constant density 
core. As the GC approaches the core, the system will need to 
rearrange itself and reach a new equilibrium state. The non- 
linear interplay between the GC and the background distri- 
bution during this rearrangement leads to a period of en- 
hanced, super-resonant friction. After ~ 1 dynamical time, 
the distribution function of the background will now be the 
correct one for the GC plus harmonic core, and dynamical 
friction will cease. Note that this rearrangement may also 
be understood in terms of the TW84 dynamical feedback 
discussed in section ITTI 

Our model does not include the back-reaction of the 
test particles on the GC, nor does it include the interaction 
between the background particles themselves. However, we 
find a good agreement between our analytic model and full 
N-body simulations, which include the above effects, in sec- 
tion|3 This suggests that our simple model does capture the 
essential physics of the problem. 

The analytic set-up is shown in Figure The in-falling 
GC at a radius, is marked by the black circle and is a 
phase angle, a, away from a given background particle at 
a radius, r^. We assume that the underlying potential is 
always harmonic (given by equation |7|l ; the GC is well ap- 
proximated by a point mass; and the background potential 
is nailed down (this is reasonable provided that Mc <^ Me„ 
where Me„ is the mass enclosed by the GC). Under these as- 
sumptions, the equation of motion for a single background, 
massless tracer, particle is given by: 
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Figure 1. A schematic diagram of the analytic set-up. The in- 
falling GC is marked by the solid black circle. The GC and particle 
orbits are marked by the grey ellipses. See text for further details. 



r -\- \ l r = r = — 



(10) 



where F_ is the specific force on the particle from the GC, G 
is the gravitational constant and Mc is the mass of the GC. 

We now search for stable solutions to eatiation llOl where 
the GC orbit is unchanged by the background. Combining 
equations |5] and IIUI gives : 



- 



(11) 



where r^ = r_^ — rp- 

From equation 1111 it is clear that stable solutions ex- 
ist where the background particles move on circular epicy- 
cles about the GC with \r^\ = const. However, more gen- 
eral solutions may be found by noting that equation 1111 
is spherically symmetric. Moving to spherical polar coor- 
dinates: — {r,9,(j)), it is straightforward to show that 
the 6 and specific angular momenta are conserved: Jg = 
r^9 = const; = const. This is to be expected given the 
symmetry of the problem. Equation II II then reduces to: 



dr 



= 



GM, ^ J| 
r 2r2 



(12) 



(13) 



where $cff is the effective potential. 

From eauation ll2l we can see that in general, the back- 
ground particles move on epicycles about the GC. These 
epicycles will not be closed, but they are quasi-periodic. Pro- 
vided the distribution function of these background particles 
is the correct equilibrium distribution for the GC plus the 
harmonic core, there will be no time averaged momentum 
exchange between the GC and the background, and there- 
fore no dynamical friction. The epicyclic orbits are stable, 
as can be readily seen by considering ^^5^. 

The existence of stable analytic solutions is a very spe- 
cial property of the harmonic potential; they exist because 
Q = const. In more general potentials, Q = 0(rp) and the 
symmetry of equation 1111 is broken. This is an important 
issue. One can imagine a thought experiment where a GC is 
held (artificially) on a fixed orbit in a general spherical po- 
tential. After a few dynamical times, it will have scattered 
the resonant background particles, reducing the torque from 
the background to zero. It is important to stress that this 
is quite different to the situation we have described in this 



paper. In the above thought experiment, a tiny perturbation 
to the GC orbit (which must in practise occur as a result of 
its self-consistent interaction with the background) will ex- 
pose the GC to an entirely new set of resonant background 
particles: dynamical friction will not cease. In our example, 
however, any perturbation to the GC orbit will not alter 
its orbital frequency at all: the resonances will remain un- 
changed. This is why our assumption, above, that the GC 
orbit is fixed is not an important one for the harmonic po- 
tential, but would be for any other potential. We test that 
this is indeed the case by relaxing the assumption of a fixed 
GC orbit in sectional 

We can use the above solution to calculate the final dis- 
tribution of background particles at equilibrium when the 
dynamical friction ceases. The key point is that the final 
distribution will move on stable epicycles about the GC. 
Firstly, this means that we can expect a density enhance- 
ment around the GC, and a depletion of particles away 
from the GC. Secondly, we can expect a large depletion in 
counter-rotating particles with respect to the centre of the 
potential. All particles, whether they move on co- or counter- 
rotating epicycles have guiding centres which co-rotate with 
the GC. For \rj < \rj, the radius of the epicyclic orbit 
is smaller than that of the GC: none of these particles can 
counter-rotate with respect to the centre of the potential. 
For |r^| > |r^|, particles on counter-rotating epicycles can 
appear to counter-rotate with respect to the centre of the 
potential. These will be a small fraction of the total particles 
which remain. We test these qualitative expectations, using 
simulations, in sectional 

A final point, which will become important later on, is 
that the orbit plane of the GC matters. Eauation llll is spher- 
ically symmetric about the GC and hides this fact. If the GC 
orbit changes (and noise within the full N-body simulations 
can cause this to happen) then the angular frequency vec- 
tor of the GC, f2, will change: the background distribution 
will no longer be in equilibrium with the GC. The system 
will have to move once again into equilibrium and this re- 
arrangement will lead to some associated dynamical friction 
on the GC. 



2.3 The Kalnajs solution 

Our analytic m ethod is a more general case of an earlier 
result found bv iKalnai j I^T^)- Kalnajs studied dynami- 
cal friction in a uniformly rotating sheet in which all par- 
ticles initially move on circular orbits. This is an equiva- 
lent problem to a GC moving on a circular orbit within a 
harmonic potential. He showed, using results from plasma 
physics, that in this case dynamical friction will vanish. Here 
we generalise this result to a GC moving on a general orbit 
within a harmonic potential. In our solution, the background 
perturbation need not lie in the plane of the GC orbit. 



3 SIMULATIONS 

In this section we compare semi-analytic and full N-body 
simulations to the analytic formulae derived in section |5| 
The simulations are labelled as in Table Q and described in 
detail in the sub-sections below. 

The analytic arguments given in section |5| suggest that 
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once a GC is introduced to a constant density background, 
the system will move towards a stable equilibrium where 
the background particles move on epicycles about the GC. 
However, this simple analytic argument cannot say anything 
about interactions between the GC and the background par- 
ticles prior to such an equilibrium being achieved; or of in- 
teractions between the background particles themselves. In 
this section we investigate this approach to equilibrium using 
numerical simulations. We use two types of simulation. The 
semi-analytic (SA) run solves equation IIUI numericallv. We 
still assume that the GC orbit is fixed, however we can study 
how the system moves from one equilibrium state (without 
the GC) to its final equilibrium with the GC. The full N- 
body (NB) run includes the interaction between the particles 
and the GC self-consistently. The GC is now free to respond 
to the background particles. This allows us to study the full 
effect of dynamical friction on the GC as the system moves 
towards its new equilibrium state. We compare results for a 
GC initially on a circular orbit and an elliptical orbit. 



3.1 The semi-analytic model: SA 

In the semi-analytic model (SA) we solved equation 1101 
with the GC orbit held fixed (the GC initial conditions 
are described in section [3.3(1 . The underlying potential was 
pure harmonic and static. We used = 4/37rGpo and 
Po = 9.93 X 10^ Mq kpc~'^. We used an isotropic, constant 
density, 3D, initial distribution of massless tracer particles. 
The equations of motion were solved using an RK4 numer- 
ical integrator ( Press et al.i ll992'l. with fixed timesteps of 
1.5 X 10"'' Gyrs. This was found to conserve energy to ma- 
chine accuracy over the whole simulation time in the limit 
Adc 0. We ran the simulations for 1 Gyr, which is ~ 10 
dynamical times for the GC at the core radius. This is the 
appropriate length of time for comparison with the full N- 
body run (see section lS^ . We tried runs with force softening 
for the GC and without. There was no significant change in 
the results for a force softening of 10 pc. The GC orbit was 
chosen to match the final stalled orbit observed in the N- 
body models. We ensured that the final position of the GC 
was identical in both models. 



3.2 The N-body model: NB 

In the full N-body (NB) model, we used the parallel multi- 
steppi ng N-body tree-code, pkdgrav2, developed by ■Stadel 
(|2001j). The potential was calculated self-consistently from 
the live particle distribution. The GC was allowed to freely 
respond to the background particles. 

We constructed stable particle halos using the tech- 
niques developed bv iKaza ntzidis et al. ( 200^. The particles 
are drawn self-consistently from a numerically calculated 
distribution function. We used a den sity distributio n that 
is described b y t he a, 7 law (iHernauist 1990 , Sah3[T993. 
iDehnenllTgg^ and IZhadll996l) : 



p{r) 



po 



(14) 



where we used po — 9.93 x 10^ kpc""^, rs — 0.91 kpc, a 
= 1.5, P — 3.0 and 7 — 0.0. Note that is the scale radius, 
not the core radius. The radius at which the log-slope of the 




0.01 



0.10 



1.00 



10.00 



r(kpc) 



Figure 2. The density distribution for the background particles 
used in the numerical simulations, see equation 1141 The dotted 
line marks the asymptotic central core where the density is con- 
stant and the potential harmonic. Inset in the plot is the distri- 
bution of orbital frequencies in the core region, plotted as \ri/ri\; 
Ti = Xp,yp,Zp (solid, dotted and dashed lines). These are equal 
and strongly peaked around a single value, showing that the core 
is indeed harmonic. 



density profile is shallower than —0.1 is rcore ~ 200 pc, which 
defines the constant density region in this model. This halo 
has a virial mass of 2.0 x 10^ Mq and the concentration pa- 
rameter is 40. A plot of the density profile is given in Figure 
121 where rcoro is marked by the vertical dotted line. Inset in 
the plot is the distribution of orbital frequencies in the core 
region, plotted as |ri/ri|; n — Xp,yp,Zp (solid, dotted and 
dashed lines). These are equal and strongly peaked around 
a single value, showing that the core is indeed harmonic (c.f. 
equation . 

The NB run, with lO'^ particles, corresponds to just 
10'^ particles within 300 pc. To achieve higher resolution, 
in the NB3 model we also used a novel 3-shell approach 
iZemp et alj |200(J) . We briefiy summ arise this approac h 
here, but defer the details and tests to lZemp et al.l i200£t) . 
The 3-shell model breaks up the mass distribution into three 
concentric spheres. The particles in each sphere are reduced 
in mass and increased in number so that central regions are 
of higher resolution. Such a model is very useful for the cur- 
rent study where we would like many particles to accurately 
sample the central harmonic core, but are not interested in 
the outer density profile which may then be less accurately 
sampled. Massive particles from the outer sphere can and 
do enter the central core in this model, but they are given 
proportionately higher force softening to prevent them from 
causing spurious hard scattering. The model produces sta- 
ble density profiles over > 20 Gyrs, very high central resolu- 
tion, and no unwanted two body effects. More detailed tests 
are given in IZemp et alj (l200dl . but for the present study 
we also explicitly verified that the single component model 
(NB) gives comparable two-body noise (see appendix IXt. 
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Simulation 


Description 


GC orbit 


Potential 


Resolution 


Mc 


SA{c,e} 


Semi-analytic 


Fixed, {(c)irc., (e)llip.} 


Fixed, harmonic 


10^ tracer 


2 X IO^Mq 


NB{c,e} 


N-body 


Live 


Live, Q,/3,7 = [1.5,3,0] 


10^ 


2 X lO^M© 


NB3{c,e} 


N-body 


Live 


Live, q:,/3,7 = [1.5,3,0*] 


10^ 3-shell 


2 X 10-^ 



Table 1. Simulation labels and parameters. From left to right the columns show the simulation label; a brief description (for more details 
see the relevant subsection in section |3J; the initial GC orbit (see section l3.3i : the background gravitational potential; the simulation 
resolution; and the mass of the GC (Mc). Parameters marked with a * are allowed to vary. In section 14.31 we measure the effect of 
changing 7 on the NB3 simulation; in section f4. 21 we measure the effect of changing Mc. 



We used a three shell model that has 10 particles for 
the innermost sphere with 300 pc radius, 10® particles for 
the shell between 0.3 and 1.1 kpc and 4 x 10® particles for 
the rest of the halo. This gives us 0(10®) particles within the 
core region. To achieve a similar number of particles within 
the central 300 pc without the 3-shell model would require 
4 X 10* particles in total. This is not yet technically feasible. 
Yet, as we show in section [3.41 and appendix 1X1 such high 
resolution is required to avoid spurious precession of the GC 
orbit plane. The advantage of the 3-shell model, given such 
limitations, is clear. The softening lengths of the particles in 
these shells were 3, 30 and 300 pc respectively. The particle 
masses were 8.9 Mq, 164.0 Mq and 757.2 Mq. Even the most 
massive particles were 100 times less massive than the GC. 
We experimented with varying the shell force softening and 
radii and found our results to be insensitive to these values. 



3.3 The GC orbit 

We used a GC mass of Mc = 2 x 10^ Mq with a force soft- 
ening of 10 pc. This gives Mc/M^n = 0.06, where Men is 
the total mass of background particles inside rcoro. For the 
N-body simulations, the GC was placed initially at a ra- 
dius of 1.069 kpc on a circular (NB3c) and elliptical (NB3e; 
Vi = 0.4vcirc) orbit. In both cases the GC orbited in the 
Xp,yp plane. For the SA simulations, the GC orbit was cho- 
sen to match that of the GC in the N-body models after it 
hit the constant density core and stalled. Its orbital phase 
was chosen such that at the end of the SA simulation (after 
1 Gyr) the GC would be in the same place as in the N-body 
simulations. 



3.4 Particle noise, resolution and convergence 

We had surprising difficulty in obtaining enough resolution 
in the N-body simulations for our results to be believable. 
The problem centred around the precession of the GC orbit 
plane. For a spherical potential (such as that studied here) 
all orbits, including that of the GC, should be planar. How- 
ever, in our initial lower resolution runs, with a resolution of 
10^ particles within 300 pc, we found that the GC orbit plane 
would precess, sometimes by as much as 20° over 10 Gyrs. 
Since we are trying to model an effect that relies critically 
on the orientation of the GC orbit plane, it is essential that 
the plane remains stable. 

In appendix ^ we use a simple analytic model of a 
2D random walk to prove that this precession is a result 
of two-body noise in the simulations. Reducing such noise 
drove us to use the three shell model discussed above. We 



show in appendix^ that the noise is not some special prop- 
erty of the three shell model, but is present in all N-body 
simulations. We found that some initial GC orientations 
showed more precession than others, for the same resolu- 
tion. This is to be expected from a random walk driven by 
two-body noise. The effect of such precession was found to 
be quite small. However, it does lead to a spurious (and 
very slow, sub-Chandrasekhar) decay of the GC orbit once 
it reaches the core. We present the results here from sim- 
ulations which showed the minimal GC plane precession. 
However, our main results are not sensitive to such selec- 
tion. Nor are our results sensitive to the use of the three 
shell model. 



4 RESULTS 

4.1 The stalling of dynamical friction in the core 

Figure |5] (straight solid line) shows the decay of the radius 
of the GC as a function of time for the NB3c and NB3e 
simulations (see Table0 . Overlaid is the prediction from the 
Chandrasekhar formula given in equation^ For this we used 
a constant InA = 5, which is the value we use throughout 
this paper. If we equate bmin with the GC force softening, 
femin = 10 pc, this gives feniax — 1.5 kpc, which is of order the 
'size' of our system. This is consistent with values found in 
other numerical s tudies of dynam ical friction on point mass 
particles (see e.g. ISpinnato et al ]|2003). 

As the cluster nears the constant density core (rcore ~ 
200 pc), it enters a phase of super-Chandrasekhar dynamical 
friction, after which dynamical friction practically ceases. 
This occurs irrespective of the initial GC orbit. This is in 
excellent qualitative agreement with analytic expectations 
from section |5| 

Figure 2] shows the distribution of particles in a slice 
about the orbit plane of the GC. The slice is defined such 
that I Jp-Jcl < li^pliil cos(6l), with 6 = 10°, where J^^ is the 
specific angular momentum of the particle and GC respec- 
tively. The left panel shows density contours for the parti- 
cle distribution (which was initially constant-density) in the 
Xp,yp plane. The right panel shows velocity histograms for 
the component of the velocity; where i;^ is the velocity 
about the Zp-axis. We do not show the Vr and vq compo- 
nents of the velocity, since they are not altered from the 
initial conditions and remain approximately Gaussian (r, 9 
and are the usual spherical polar coordinates). In the top 
panels, the solid lines show the slice just before the GC hits 
the core in the NB3c simulation (at time t = 5 Gyrs); the 
dotted contours show the SAc simulation at t = 0. The mid- 
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Figure 3. The decay of the radius of the GC as a function of 
time for a GC on a circular (straight solid line) and elliptical 
(oscillating solid line) orbit. Overlaid (dotted lines) are the pre- 
dictions from the Chandrasekhar formula given in equation 
using In A = 5. Notice that for the first few Gyrs the agree- 
ment with equation is excellent. As the cluster nears the con- 
stant density core (rcoro ~ 200 pc), it enters a phase of super- 
Chandrasekhar dynamical friction, after which dynamical friction 
practically ceases. 



die panels show similar results for the NB3c and SAc sim- 
ulations at times t — 8 Gyrs and t — 1 Gyrs, respsectively. 
The bottom panels show the NB3e and SAe simulations at 
times t — 4 Gyrs and t = 1 Gyrs, respectively. We analyse 
within this slice to highlight the changes in density caused 
by the GC. Outside of the slice, background particles still 
move on epicycles about the GC, but their projected posi- 
tions onto the Xp, yp plane make it difficult to see the density 
enhancement about the GC. 

Notice that the velocity histograms for (right panel, 
top) are double-peaked. This is because we have taken a thin 
slice in the orbit-plane of the cluster. The only particles that 
have v,^ = in this plane are those on pure radial orbits, 
which is a very small number of particles in the isotropic 
initial conditions. 

In the top panels of Figure |1] the particles are close to 
their initial configuration. There has been some depletion in 
density at the centre and the on-set of some substructure, 
but the velocity histograms show that the velocity distribu- 
tion of the background particles is still isotropic. 

The middle and bottom panels in Figure 01 show 
the distribution of particles in the slice after the super- 
Chandrasekhar friction has ended, and the GC has settled 
into equilibrium in the core. Notice the good agreement with 
the semi-analytic (SAc,e) simulations for both the density 
and velocity distribution in the slice, irrespective of the ini- 
tial GC orbit. As expected from the arguments given in sec- 
tion 121 the number of counter-rotating particles has been 
significantly depleted. 

The density distribution in the slice is peaked just be- 



hind the cluster; it has a tail which is longer for the full 
N-body run. This is likely due to particle-particle scattering 
which prevents high density regions from forming. Such a 
tail should lead to some dynamical friction on the GC from 
the background. We estimated the strength of this effect for 
the SAc model. To do this, we summed the force from all 
of the background particles on the GC, assuming that their 
total mass was Men- What really matters is the time aver- 
aged force on the GC. This must be small, since little or no 
dynamical friction is observed after the GC reaches the core. 
However, even the total force at an instant is always smaller 
than the dynamical friction force, computed from equation 

m 

Notice from equation 1121 that we could construct any 
final density distribution using an appropriate combination 
of epicyclic orbits about the GC. In practise, however, the 
final density distribution is set by the initial configuration 
of background particles within the core. The transformation 
of this initial distribution by the arrival of the GC, must be 
determined numerically. The SA model is essential in this 
respect. 

The keen observer will notice that the enhanced friction 
appears to set in rather near the region where the resolution 
in the 3-shell model increases (recall that the high resolu- 
tion inner shell starts at 300 pc) . This is almost certainly a 
coincidence. We performed two tests to check this. Firstly, 
an explicit test by starting a GC sinking inside the high 
resolution shell. Once again, we observed enhanced friction 
followed by stalling. Secondly, we performed a test-run start- 
ing the GC outside the 2nd shell. As it moved through the 
shell transition at 1.1 kpc, no detectable effect was observed. 

4.2 The effect of varying Mc 

Figure |^ shows the decay rate of the GC as a function 
of the GC mass, Mc (solid lines); Mc is marked in solar 
masses. Overlaid are analytic predictions from equation Q 
using In A = 5, as previously (dotted lines). For these sim- 
ulations, we re-ran simulation NB3c but using 10 times the 
GC mass {Mc = 2 x IO'^Mq), and half of the GC mass 
{Ale = 1 X 10"" M0). For reasons of computational expense, 
we ran these new runs at a lower resolution with 10® par- 
ticles in each shell. We could not investigate smaller GC 
masses than Mc = W^Mq, since then Mc approaches the 
mass of the heaviest particle and two-body effects dominate 
over dynamical friction. 

Notice that the lower resolution runs are noisier and 
decay faster once the GC hits the core. This decay is due 
to the precession (due to numerical noise) of the GC orbit 
plane, discussed in section and appendix^ and is much 
smaller in the higher resolution runs. We explicitly checked 
that this is indeed the case using lower resolution runs of 
NB3c. 

In all runs, the GC shows a reduced friction at the core 
region. Notice that the point at which the GC departs from 
Chandrasekhar like friction appears to be a weak function 
of the GC mass. This is to be expected: a more massive GC 
will more rapidly scatter the background particles and stall 
more quickly once it reaches the core region. However, it is 
tempting to suggest a simpler explanation: that equation Q 
is failing simply because Me„ = r/Mc; where Men is the final 
mass enclosed and ri is some constant of order unity. This 
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Figure 4. The distribution of particles in the Xp,yp plane for the N-body (solid lines) and semi-analytic (dotted lines) simulations. 
The left panels show density contours for the particle distribution (that was initially constant-density) in the Xp,yp plane. The right 
panels show velocity histograms for the component of the velocity; where is the velocity about the Zp-axis. The solid vertical line 
marks the GC velocity about the Zp-axis. The top two panels show the GO circular orbit simulations just before the GC experiences 

super-Chaiidrasckliar dynamical friction. This corresponds to t = OGyrs for the SAc simulation, and t = 5Gyrs for NB3c. Notice that, 
for the NB3c sinmlation (solid lines), the background particles are nearly unchanged from their initial distribution. The middle two 
panels show the GC on a circular orbit after the super- Chandrasckhar friction has finished and the GC has settled into a steady-state 
in the harmonic core. This corresponds to t = 1 Gyrs for the SAc simulation, and t = 8 Gyrs for NB3c. The bottom two panels show the 
GC on an elliptical orbit {vi = 0.4ticirc) after it has reached the harmonic core. This corresponds to t = 1 Gyrs for the SAe simulation, 
and t = 4 Gyrs for NBSe. In all cases the final position of the GC in the SAc,e and NB3c,e simulations is identical and marked by the 
solid circle. Note that there is no GC marked in the top panels since, in the NB3c simulation, the GC lies outside of the plot area at this 
time. 
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Figure 5. The decay rate of the GC as a function of Mc (soUd 
lines); Mc is marked in solar masses. Overlaid are analytic predic- 
tions from equation^using In A = 5, as previously (dotted lines). 
For these simulations, we re-ran simulation NB3c (see Table 
but using 10 times the GC mass (Mc = 2 X lO'^M©), and half 
of the GC mass (Mc = 1 X 10^ M©). Me„/Mc at the point of 
the onset of the stalling behaviour is, in order of increasing Mc'. 
[3,5,2], 



Figure 6. The decay rate of the GC as a function of the central 
log-slope of the background density distribution, 7 (solid lines); 
7 is marked on the plot. Overlaid are analytic predictions from 
eouationin using In A = [8, 7, 3.5], in order of increasing 7 (dotted 
lines). The crosses mark the radii at which the final density profile 
has a central log-slope shallower than -0.1. M^n/Mc at the point 
of the onset of the stalling behaviour is, in order of increasing 7: 
[4,1]. 



is perhaps worrisome given the mass ratios, M^n/Mc, at the 
point of the onset of the stalling behaviour. These are, in 
order of increasing Mc- [3,5,2]. However, we believe that the 
situation is not this simple for the following reasons: (i) if the 
stalling were a result only of Mc — Men, then it would not 
be a special property of constant density cores. We show in 
sect ion 1^31 below, that the stalling behaviour does not occur 
for steeper density profiles, whatever the enclosed mass, (ii) 
The model we present in section |5| provides a good fit to 
the final density and velocity distribution of the background 
particles in the core suggesting that we have captured the 
correct physical explanation. 

4.3 The effect of varying 7 

Figure|^shows the decay rate of the GC as a function of the 
central log-slope of the background density distribution, 7 
(solid lines); 7 is marked on the plot. For these simulations, 
we re-ran simulation NB3c but using 7 = [0.1,0.5]. Also 
shown are re s ults f or a simulation with 7 = 1 taken from 
iGoerdt et alJ l)200(fl . Overlaid are analytic predictions from 
equation using InA = [8,7,3.5], in order of increasing 7 
(dotted lines). InA is different for each of these simulations, 
reflecting the change in the underlying density distribution; 
similar results have been found elsewhere in the literature 
(see e.g. l.Tust fc PenarrubialEoOSh . All simulations were high 
resolution (~ 10^ particles per shell), but since we are inter- 
ested in the core stalling properties of the GC, we started the 
7 = [0.1,0.5] simulations at ~ 400 pc, rather than ~ Ikpc 
as previously. 

The key point is that the 7 = 1 model is well-fit by the 



Chandrasekhar form over the entire simulation time. This 
is despite the fact that Men — Mc at ~ 0.1 kpc for this 
run. This suggests that the core stalling behaviour is a spe- 
cial property of the harmonic core and not to do with the 
enclosed mass. However, the 7 = [0.1,0.5] runs both show 
stalling behaviour despite not having a central core. This 
occurs because the GC itself creates a small core as it falls 
in and heats the background particle distribution. For initial 
density distributions steeper than 7 = 0.5 this no longer oc- 
curs. In this case, the density profile does become shallower 
as a result of heating, but the heating is not sufficient to 
form a core in the centre before the GC falls all of the way 
in. The crosses on Figure |S| mark the radii at which the final 
density profile has a central log-slope shallower thant -0.1. 
Recall that this is the same definition we used to define rcoro 
earlier. 



5 CONCLUSIONS 

Using analytic calculations and N-body simulations we have 
shown that in constant density harmonic cores, sinking 
satellites undergo an initial phase of very rapid (super- 
Chandrasekhar) dynamical friction, after which they expe- 
rience no dynamical friction at all. This occurs because, for 
the special case of harmonic potentials, there are stable so- 
lutions where the background particles move on epicycles 
about the in-falling satellite. The system moves rapidly into 
this stable configuration. In doing so, the satellite experi- 
ences a brief moment of enhanced friction. Once in equilib- 
rium, there is no net momentum transfer between the back- 
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ground particles and the satellite and friction ceases. For 
density profiles with a central power law profile, p oc r'" , 
the infalling satellite heats the background and causes a to 
decrease. For a < 0.5 initially, the satellite generates a small 
central constant density core and stalls as in the a — case. 

Our results concerning dynamical friction stalling in 
constant density cores are of broad astrophysical interest. 
Recent observational work suggests that galaxies may have 
central dark matter density cores, rather than the den- 
sity cusps predicted by numerical simulations. Galactic bars 
orbiting in such potentials will experience very weak dynam- 
ical friction and can be very long-lived (in fact central den- 
sity distributions do not need to be pure harmonic to see this 
effe ct, low-density will also lead to very little friction - see 
e.g.toebattista fc Sellwoodll99^ and lDebattista fc SellwoodI 
l200ol) . Satellites falling into such galax;ies will stall at the 
core radius and never make it to th e centre. This poin t 
was investigated in a companion paper (iGoerdt et al .120061) . 
where we suggested that a constant density core could solve 
the 'timing problem' for the GCs in the Fornax dwarf galaxy. 
Finally, recent work on merging black holes suggests that 
they can form a central cons tant density core in the b ack- 
ground distribution (see e.g. iMilosavlievic et alJl20o3 and 
iRavindranath et alj|2003) prior to forming a hard binary. If 
true, our results suggest that this could further exacerbate 
the well-known problem of getting the binaries to coalesce. 
Their rate of hardening will stall, even before the majority 
of stars and dark matter have been ejected from the core, 
if the background distribution is close to constant density. 
This may point towards gas playing a more important role 
in bringing superm assive bla ck holes together at the centre 
of galaxies (see e.g. iGould fc Rixii200g") . 
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APPENDIX A: TWO BODY NOISE AND 
PRECESSION OF THE GC ORBIT PLANE 

In this appendix we present a simple analytic model for the 
precession of the GC orbit plane due to particle noise and 
compare this with the simulations. We show that even quite 
small particle noise can lead to significant plane precession 
over ~ 100 dynamical times. 

Under the assumption of linear background particle 
trajectories, it is straightforward to show that an interac- 
tion with one background particle will produce a veloc- 
ity kick perpendicular to the GC's orbit plane given by 
jBinnev fc Tremainelll987[) : 
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Figure Al. Evolution of the GC orbit plane in angle, 6, over the 
simulation time. The straight solid lines are for an analjrtic model 
that assumes a 2D random walk. Results are shown for increasing 
particle number, A'^. Over-plotted are results from three typical 
N-body simulations: NB, NB3c and NB3c' - see Table [Tl and this 
appendix for details; NB3c' is identical to NBSc, except that the 
GC initial orbit plane is different. 



where m is the mass of the background particle, Mc is the 
mass of the GC and b is the impact parameter (initial per- 
pendicular separation) of the encounter. Such a kick occurs 
over ^ a dynamical time. 

Summing over all such encounters (all impact parame- 
ters) then gives the mean total velocity kick to the GC in 
~ a dynamical time. We sum over Sv'^ to give the r.m.s. 
change; Svz is of random sign and will sum to zero; 
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where v is now ~ the velocity of the GC; A — bmax/bmin 
is the term inside Coulomb logarithm that also appears in 
equation A'^ is the number of background particles inside 
~ &max (the GC orbit is assumed to lie in the x-y plane); 
and In A' = l/2[...] is defined by equation IA2I 

Notice that in the limit of large impact parameters, 
femax > &min > GMc/v"^ si^max > Xmin > 1, and equation 
IA2l reduces to the more famiUar form; Av'^Jv^ = {8/N) In A. 
It is then independent of the GC mass. 

In one orbit, the GC will move a mean z distance, Az ~ 
Aiiztdyn, where tdyn is the orbit time. The mean change in 
angle over one orbit, A6, of the vector normal to the GC 
orbit plane is then given by; 
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where we have assumed that the GC moves on a circular 
orbit of radius, r. 

Any dependence on the underlying potential completely 
factors out in equation IA4I and Ad depends only on the 
number of particles, A'', and very weakly on Mc/m. 

The orbit plane can be tilted due to such scattering 
noise from the background distribution in two independent 
directions. Since the potential is spherical, there is no restor- 
ing force and once the plane has tilted, the probability it will 
tilt again is independent of its past history. Thus we may 
model the accumulated precession of the orbit plane by a 
2D random walk. This gives: 

e = M.f^ (A5) 

The orbit time at r = rcorc for our model is given 
by tdyn = 2tt^J^^^ = 0.15 Gyrs. In Figure |^ we plot 
the mean orbit plane precession predicted by this random 
walk model, as a function of simulation time, t/t^yn- We 
use 6max = 1.5 kpc and bmin ~ 10 pc, which gives InA = 5, 
as in section 0] In section 31 we typically ran our N-body 
models for 10 Gyrs which corresponds to ~ 100 dynamical 
times. The straight solid lines show the effect of increasing 
the particle number, N. Notice that extremely high reso- 
lution is required to keep plane precession to a minimum 
over our simulation time: even with 10^ particles we can ex- 
pect a mean precession over the whole simulation of ~ 7°. 
Over-plotted are results from the NB, NB3c and NB3c' sim- 
ulations (see Table 0. Recall that the NB model was a sin- 
gle shell model with 10^ particles in total, with 10^ within 
300 pc. The NB3 simulations were three-shell models with 
10^ particles within 300 pc. NB3c' is a simulation which is 
identical to NB3c but with a different GC initial orbit plane. 
Notice that in all cases the plane precesses; it is not some nu- 
merical error introduced by the three shell model. The NB3 
simulations show a smaller precession than the NB simu- 
lation as is expected given their higher effective resolution. 
Finally, notice that changing the initial GC orbit plane can 
alter the total precession quite dramatically (compare the 
NB3c and NB3c' simulations). This is to be expected given 
the random walk model, above. 

The total particle number, N, in equation IA4I is a 
slightly ill-defined quantity and so should not be equated 
exactly with the number of particles in the simulation (par- 
ticularly for the 3-shell models). However, it is encouraging 
that our simple random walk model produces the correct 
mean slope for the plane precession and the correct scaling 
with particle number. It is clear that the 3-shell model has 
an advantage over the single shell model: it samples the core 
region with 1000 times the resolution of the single shell and 
shows much smaller two-body noise. 

Throughout this paper, we present simulations which 
minimise the evolution of the orbit plane. It is important 
to note, however, that all of our simulations show the same 
central result: a period of super-Chandrasekhar friction, fol- 
lowed by stalling at the constant density core. 
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